\documentclass{article}%article类型文章
\usepackage{ctex}%支持中文编译的宏包，编译一定要选择XeLatex啊！！
\usepackage[utf8]{inputenc}

\title{柱壳的振动特性分析}%标题
\author{BY2104502 杜晨鸿}%作者
\date{2021-12}%日期
 
%——————以下宏包功能不一定用上，请自行学习——————————
% \usepackage{natbib}
\usepackage[table,xcdraw]{xcolor}%提供更多的色彩
\usepackage{extarrows} %长等号扩展
\usepackage{booktabs}%三线表格
\usepackage{float}%浮动体宏包
\usepackage{amsmath}%扩充数学富符号
\usepackage{multirow}
% \usepackage{cite}
\newcommand{\upcite}[1]{\textsuperscript{\cite{#1}}}
%设定了引用命令\upcite来右上角引用
\usepackage{graphicx}%引入图片宏包
\usepackage{subfigure}%小图片
\usepackage{listings}  %代码块当然数学作业也用不到
\usepackage{makecell}

%——————————————————————————————————————————————————
\newcommand{\res}[1]{\mathrm{Res}\left[{}#1 \right]}
\newcommand{\ykh}[1]{\left({}#1\right)}%定义命令\ykh{}来输出圆括号，如果直接输出()会导致固定的括号大小与公式不匹配
\usepackage{geometry}%设定版面格式
\geometry{
  a4paper,
  top=25.4mm, bottom=25.4mm,
  headheight=2.17cm,
  headsep=4mm,
  footskip=8mm
}%设置标准页边距
\usepackage[justification=centering]{caption}%图片标题居中
%不熟悉latex的数学公式代码可以去www.latexlive.com 生成。
% 设置题注格式
\renewcommand {\thefigure} {\arabic{figure}}


\begin{document}%正文部分
\captionsetup{labelformat=default,labelsep=space}  % 去除题注冒号
\maketitle%显示正文之前的属性如作者
\tableofcontents%生成目录
\newpage%新的一页

\section{选定几何及材料参数}

柱壳的几何及材料参数见表\ref{materal}。

\begin{table}
  \centering
  \caption{几何及材料参数}
  % \begin{center}
  \begin{tabular}{l|l|l}
    \Xhline{1pt}
    参数     & 符号   & 值                         \\\hline
    壳厚     & $h$    & 1mm                        \\ \hline
    长度     & $L$    & 100mm                      \\\hline
    半径     & $R$    & 50mm                       \\\hline
    密度     & $\rho$ & $\rm 8\times10^{-9}t/mm^3$ \\ \hline
    杨氏模量 & $E$    & 200000MPa                  \\ \hline
    剪切模量 & $G$    & 76923MPa                   \\\hline
    泊松比   & $\nu$  & 0.3                        \\\hline
    \Xhline{1pt}
  \end{tabular}
  \label{materal}
\end{table}


\section{有限单元法求解}

\subsection{有限元模型}
本文所采用的有限元模型如图\ref{modal-fem}所示。采用ANSYS中的8节点四边形厚壳单元SHELL281。
共有单元31400个，节点94829个。对左侧一圈节点施加固定约束，采用分块兰索法进行模态分析。
\begin{figure}
  \centering
  \includegraphics[height=150pt]{modal-fem.png}
  \caption{有限元模型}
  \label{modal-fem}
\end{figure}


\subsection{计算结果}
有限单元法算得的模态中，每个共振频率都对应两个频率相同，振型相差一定角度(周向)的振型。
本文仅记录前三节互异模态(即频率和振型均不相同)，如图\ref{level1-fem}-图\ref{level3-fem}所示，前三阶振型分别为6节径、8节径、4节径。其频率见表\ref{freq-fem}。\\
\begin{figure}
  \centering
  \includegraphics[height=150pt]{./cyc000.jpg}
  \includegraphics[height=150pt]{./cyc001.jpg}
  \caption{有限元法一阶模态}
  \label{level1-fem}
\end{figure}

\begin{figure}
  \centering
  \includegraphics[height=150pt]{./cyc003.jpg}
  \includegraphics[height=150pt]{./cyc002.jpg}
  \caption{有限元法二阶模态}
  \label{level2-fem}
\end{figure}

\begin{figure}
  \centering
  \includegraphics[height=150pt]{./cyc004.jpg}
  \includegraphics[height=150pt]{./cyc005.jpg}
  \caption{有限元法三阶模态}
  \label{level3-fem}
\end{figure}

\begin{table}
  \centering
  \caption{有限元解的前三阶互异模态频率}
  % \begin{center}
  \begin{tabular}{l|l}
    \Xhline{1pt}
    阶数 & 频率(Hz) \\\hline
    1    & 1491.7   \\\hline
    2    & 1646.4   \\\hline
    3    & 2324.5   \\\hline
    \Xhline{1pt}
  \end{tabular}
  \label{freq-fem}
\end{table}

\section{假设模态法求解}

\subsection{动能及势能}
任意连续体的动能可以表示为其中所有质点动能的积分，即\upcite{lee2015free}
\begin{equation}
  T=\frac{1}{2}\rho\int_0^L\int_0^{\rm 2\pi}\int_{-h/2}^{h/2}
  (\dot{\xi}^2+\dot{\eta}^2+\dot{\omega}^2)R{\rm d}z{\rm d}\varphi{\rm d}x
\end{equation}
其中$x$、$\varphi$、$z$分别表示轴向、周向、径向坐标，其中$\varphi$的量纲为弧度，其余为长度。
$\xi$、$\eta$、$\omega$分别表示轴向、周向、径向位移，量纲均为长度。
对于薄壳，可假径向各质点具有相同的速度，即位移仅为$x$、$\varphi$的函数，与$z$无关，
于是可对其内层积分进行计算，将上式化简为
\begin{equation}
  T=\frac{hR\rho}{2}\int_0^L\int_0^{\rm 2\pi}
  (\dot{\xi}^2+\dot{\eta}^2+\dot{\omega}^2){\rm d}\varphi{\rm d}x
\end{equation}

任意线弹性连续体的变能可以表示为
\begin{equation}
  U=\frac{R}{2}\int_0^L\int_0^{\rm 2\pi}\int_{-h/2}^{h/2}
  (\sigma_x\epsilon_x+\sigma_{\varphi}\epsilon_{\varphi}+\sigma_z\epsilon_z
  +\tau_{xz}\gamma_{xz}+\tau_{z\varphi}\gamma_{z\varphi}+\tau_{x\varphi}\gamma_{x\varphi})
  {\rm d}z{\rm d}\varphi{\rm d}x
\end{equation}
我们将其假设为平面应力问题，忽略径向的正应力及切应力，可得\upcite{lee2015free}
\begin{equation}
  U=\frac{R}{2}\int_0^L\int_0^{\rm 2\pi}\int_{-h/2}^{h/2}
  (\sigma_x\epsilon_x+\sigma_{\varphi}\epsilon_{\varphi}
  +\tau_{x\varphi}\gamma_{x\varphi})
  {\rm d}z{\rm d}\varphi{\rm d}x
  \label{stress-energy-plane}
\end{equation}
平面应力问题的物理方程为
\begin{align}
  \sigma_x        & =\frac{E}{1-\nu^2}(\epsilon_x+\nu\epsilon_\varphi) \\
  \sigma_\varphi  & =\frac{E}{1-\nu^2}(\epsilon_\varphi+\nu\epsilon_x) \\
  \tau_{x\varphi} & =\frac{E}{1+\nu}\gamma_{x\varphi})
\end{align}
将其代如式\ref{stress-energy-plane}可得：
\begin{equation}
  U=\frac{ER}{2(1-\nu^2)}\int_0^L\int_0^{\rm 2\pi}\int_{-h/2}^{h/2}
  (\epsilon_x^2+\epsilon_\varphi^2+2\nu\epsilon_x\epsilon_\varphi
  +\frac{1-\nu}{2}\gamma_{x\varphi}^2)
  {\rm d}z{\rm d}\varphi{\rm d}x
  \label{stress-energy-plane2}
\end{equation}

根据Donnell–Mushtari理论\upcite{lee2015free}，应变可以表示为
\begin{align}
  \epsilon_x       & =\frac{\partial \xi}{\partial x}-z\frac{\partial^2\omega}{\partial x^2}                                                                           \\
  \epsilon_\varphi & =\frac{1}{R}\frac{\partial \eta}{\partial \varphi}+\frac{\omega}{R}-\frac{z}{R^2}\frac{\partial^2\omega}{\partial \varphi^2}                      \\
  \tau_{x\varphi}  & =\frac{\partial\eta}{\partial x}+\frac{1}{R}\frac{\partial \xi}{\partial \varphi}-\frac{2z}{R}\frac{\partial^2\omega}{\partial x\partial \varphi}
\end{align}
将上式代入\ref{stress-energy-plane2}即可求出系统的应变能。

\subsection{刚度矩阵及质量矩阵}
根据拉格朗日法，系统自由振动的动力学方程可写为
\begin{equation}
  \frac{\rm d}{{\rm d}t}\frac{\partial T}{\partial \dot{\psi}_i}+\frac{\partial U}{\partial \psi_i}=0
\end{equation}
将$i$取遍全部广义坐标即可得到系统的动力学方程。
考虑到线性系统的动力学方程可以表达为矩阵形式，且$\frac{\rm d}{{\rm d}t}\frac{\partial T}{\partial \dot{\psi}_i}$对应惯性项；
$\frac{\partial U}{\partial \psi_i}$对应刚度项，因此必然有
\begin{equation}
  \frac{\rm d}{{\rm d}t}\frac{\partial T}{\partial \dot{\psi}_i}=\sum_{j=1}^{n}m_{ij}\ddot{\psi_j}
\end{equation}
其中,$m_{ij}$表示质量矩阵第$i$行第$j$列的元素。

同理，对于刚度矩阵也有
\begin{equation}
  \frac{\partial U}{\partial \psi_i}=\sum_{j=1}^{n}k_{ij}\psi_j
\end{equation}
其中$k_{ij}$表示刚度矩阵第$i$行第$j$列的元素。

\subsection{假设振型及边界条件}
由于位移边界条件为固定约束，固柱壳端面的位移及转角均为0。
易知$\frac{\partial \omega}{\partial x}$可表示以$\phi$为轴的旋转；
$\frac{\partial \xi}{\partial x}$可表示以$x$为轴的旋转。因此在固定端($x=0$)，应有
\begin{align}
  \xi                                & =0 \\
  \eta                               & =0 \\
  \omega                             & =0 \\
  \frac{\partial \omega}{\partial x} & =0 \\
  \frac{\partial \xi}{\partial x}    & =0
\end{align}
由于$\xi$分量只要求函数值为0，因此采用一端固支杆的振型作为其$x$方向振型函数，即
\begin{align}
  \sin( & \frac{{\rm \pi} x}{2L})  \\
  \sin( & \frac{3{\rm \pi} x}{2L})
\end{align}
由于$\eta$及$\omega$分量的函数值及一阶导数均要求为0，故选用振型函数为
\begin{align}
  x^2 \\
  x^3 \\
  x^4
\end{align}


在$\varphi$方向上，根据有限单元法的结果，柱壳呈节径型振动。由于各阶振型在周向上呈现周期性，因此应使用三角函数作为假设模态。

文献\cite{lee2015free}中推荐使用$\cos n \varphi$表示$\xi$及$\omega$分量的$2n$节径振动；
使用$\sin n \varphi$表示$\eta$分量的$2n$节径振动。本文中对有限单元法的计算结果进行检查，
发现在各阶振型中，$\xi$及$\omega$分量振动的波峰刚好与$\eta$分量振动的波谷重合，
因此文献中给出的周相位移模式是合理的。
从结构的对称性角度分析，对于子午面，$\xi$及$\omega$是对称分量；而$\eta$是反对称分量，所以这种模态假设方式也是有依据的。

有限单元解中，前三阶振型的节径数分别为6、8、4；
周期分别为$\frac{\rm 2\pi}{3}$、$\frac{\rm 2\pi}{4}$、$\frac{\rm 2\pi}{2}$。
因此，本文采用$\cos(2\varphi)$、$\cos(3\varphi)$、$\cos(4\varphi)$作为$\xi$及$\omega$分量在$\varphi$方向的振型函数；
采用$\sin(2\varphi)$、$\sin(3\varphi)$、$\sin(4\varphi)$作为$\eta$分量在$\varphi$方向的振型函数。

在实际计算中发现，不同节径数的模态是互不耦合的，即在任意一阶模态中，节径数不同的振型函数只有一个非零。
因此，本文对$\omega$分量选用以下振型函数
\begin{align}
  \cos(n\varphi)x^m
\end{align}
其中$n=2,3,4$，$m=2,3,4$。
$\eta$分量选用以下振型函数
\begin{equation}
  \sin(n\varphi)x^m
\end{equation}
其中$n=2,3,4$，$m=2,3,4$。
对$\xi$分量选用以下振型函数
\begin{equation}
  \cos(n\varphi)\sin( \frac{m{\rm \pi} x}{2L})
\end{equation}
其中$n=2,3,4$，$m=1,2$。

将$n$依次取3、4、2以计算前三阶振型。

\subsection{计算结果}
前三阶互异振型的频率见表\ref{freq-mod-ana}；其振型见图\ref{level1-ana}-图\ref{level3-ana}。
\begin{table}
  \centering
  \caption{假设模态解的前三阶互异模态频率}
  % \begin{center}
  \begin{tabular}{l|l|l}
    \Xhline{1pt}
    阶数 & 频率(Hz) & 相对有限元解的误差(\%) \\\hline
    1    & 1567.8   & 5.1                    \\\hline
    2    & 1746.6   & 6.1                    \\\hline
    3    & 2403.9   & 3.4                    \\\hline
    \Xhline{1pt}
  \end{tabular}
  \label{freq-mod-ana}
\end{table}

\begin{figure}
  \centering
  \includegraphics[height=150pt]{./mod1-ana.png}
  \caption{假设模态法一阶模态}
  \label{level1-ana}
\end{figure}

\begin{figure}
  \centering
  \includegraphics[height=150pt]{./mod2-ana.png}
  \caption{假设模态法二阶模态}
  \label{level2-ana}
\end{figure}

\begin{figure}
  \centering
  \includegraphics[height=150pt]{./mod3-ana.png}
  \caption{假设模态法三阶模态}
  \label{level3-ana}
\end{figure}
\section{总结}
通过对比假设模态法和有限单元法的计算结果，可以看出，假设模态法对频率的估计与有限单元法比较接近，且误差随振动节径数增加而增大；二者的振型一致。
振动频率的估计值总是偏大，这可能是因为假设模态法强制振型为有限个给定的函数，因此模型受到了额外的“约束”，导致振动频率增加。

\bibliographystyle{gbt7714-2005}
\bibliography{ref}

\end{document}